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Abstract. 

We describe here the extension of the density matrix renormalization group algorithm to 
the case where the Hamiltonian has a non-Abelian global symmetry group. The block states 
transform as irreducible representations of the non-Abelian group. Since the representations are 
multi-dimensional, a single block state in the new representation corresponds to multiple states 
of the original density matrix renormalization group basis. We demonstrate the usefulness of 
the construction via the the one- dimensional Hubbard model as the symmetry group is enlarged 
from (7(1) x 17(1), up to SU{2) x 577(2). 



In past years, the density matrix renormalization group (DMRG) method [1] has been 
extensively used to study one and two dimensional strongly correlated electron systems [2]. 
This method became very popular when it was realized that it enabled a level of numerical 
accuracy for one dimensional systems that was not possible using other methods [3] . 

One major drawback of DMRG is that calculations are performed in a subspace of purely 
Abelian symmetries, such as the C7(l) symmetries of total particle number and the z com- 
ponent of the total spin. Thus one can only obtain a few states in different total particle 
number and z component of total spin sectors [4] . For models where ferromagnetism emerges 
the situation worsens, that is, to determine magnetization, a combination of methods must be 
employed which will artificially raise the energy of the higher spin state [5] within the chosen 
z component total spin sector. 

In recognizing the imperative need, to introduce a DMRG method which has a total spin 
quantum number naturally implemented, a number of unsuccessful attempts were previously 
made (e.g., for the spin 1 Heisenberg model [6] and t-t'-U model [7,8]). The most successful 
previous work on the application of non-trivial symmetries is the IRF-DMRG method intro- 
duced by Sierra and Nishino [9], whereby the vertex hamiltonian is first transformed into an 
interaction round a face hamiltonian [11], and then a variant of DMRG is applied to the IRF 
model. The IRF model can be chosen such that it explicitly factors out the global symmetry 
group. This technique has been successfully applied to the spin 1/2 Heisenberg chain and 
the XXZ chain with quantum group symmetry SU q (2) [9] and later, the spin 1 and spin 2 
Heisenberg chains [10]. However, the IRF-DMRG method is complicated by the necessity 
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to calculate the IRF weights for each interaction term in the Hamiltonian. The number of 
non-trivial IRF weights increases rather quickly as the magnitude of the spins in the system 
is increased and for a lager global symmetry group. Thus as far as we know, the IRF-DMRG 
has not yet been applied to any more complex models, such as a fermionic system. In the 
present work we show that non-Abelian symmetries can be naturally accommodated into 
DMRG without the need for a vertex-IRF transformation. In this form, the starting point of 
a calculation is the matrix elements of the single site operators which are relatively simple to 
calculate; the number of such elements varies inversely with the dimension of the irreducible 
representations of the global symmetry group, and thus is reduced for a larger global symmetry 
group. For example, all single site operators of a spin chain are represented as 1 x 1 matrices, 
independent of the magnitude of the actual spins. For Hubbard- type models, the number of 
single site basis states is reduced from four to two, corresponding to the spinon and holon 
states [19]. 

We do not attempt here to give a complete description of the DMRG algorithm, instead 
we refer the reader to the original description by White [1] and more recent reviews [3], and 
concentrate on the essential elements of the algorithm that require modification when using 
non-Abelian symmetries. These are the construction of tensor product basis and operators 
(whether it is through adding a single site to a block, or joining blocks to construct a su- 
pcrblock) , and the truncation of block states via the reduced density matrix. 

We introduce the method by way of the Lie group SU(2). This symmetry is readily 
applicable to all quantum spin systems that can be written in a form that does not break 
rotational symmetry. In principle, it is not difficult to calculate eigenstates of SU(2) for a 
finite system by using the Clebsch-Gordan transformation [12], especially in DMRG where 
the system is built one or two lattice sites at a time. In this construction, the tensor product 
of two basis vectors, labelled here by subscripts 1 and 2, is 

\jm(jihaia 2 )) = ^ C^ 2m |jimi(ai)) |j 2 m 2 (a 2 )) , (1) 

mi ,7712 

where C^^ 2m is the Clebsch-Gordan coefficient. Here we use the notation j is the total 
spin quantum number, S 2 \j) — j(j + l)|j), m is the projection of the spin onto the z— axis 
and (a) is an index that encapsulates the additional labels used in DMRG (ie, to label the 
a'th basis state of the given quantum numbers). Bracketed labels are not associated with a 
quantum number. Constructing basis states in this way in DMRG suffers from two problems. 
Applying this transformation involves two summations for each operator matrix element. This 
impacts severely on the computational effort required to construct the block, and especially 
the superblock, operators. Secondly, the direct application of the usual DMRG reduced 
density matrix to a wavefunction constructed from some (jm) subspace of Eq. (1) does not 
commute with the SU(2) generators. Indeed, the wavefunction of an ST/ (2) invariant system 
represents the same physical state independent of the z-component of the spin, so the density 
matrix of the full system is 

P = b' m ( a '))^ m (a')^im(a)0' m ( a )l i ( 2 ) 

m.a' ,a 

where the wavefunction \^ m ) — J2 a ipj m (a) \jm(<x)) is an eigenstate of total spin, ie S + |<I' m ) = 
\J (j — m)(j + m + l)|W m+ i). Using Eq. (1) and tracing over the right basis, the SU(2) 
invariant reduced density matrix for the left block can be constructed, 

mi("i); .727712(0:2) 'rj 1 m 1 (a 1 ); j 2 m 2 (a 2 ) ' 

(3) 

"H ,J2 ,777-2,0:2 
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which acts identically on each mi component of the basis. This equation can also be directly 
by adding an additional constraint to the reduced density matrix, to force each eigenstate of 
the reduced density matrix to be an eigenstate of S 2 . For j = this reduces to the usual 
DMRG density matrix. This is seen in conventional DMRG by the well known (2ji + l)-fold 
degeneracies in the reduced density matrix eigenvalues. 

Despite the additional overhead of the Clebsch-Gordan transformation, this construction 
of SU(2) invariant DMRG works well for small values of j, and is described further in [15]. 
However, further improvements are possible. The projection quantum number m can be 
completely eliminated using the Wigner-Eckart theorem, 



(j'm'{a')\T J M \jm{a)) = C^ m ,(j'{a')\\T J \\j{ a )) , 



(4) 



for the M'th component of an operator T J transforming as a rank 2 J+ 1 tensor. The quantity 
(j'(o ! ')||T' J |b( Q; )) is the reduced matrix element [12] and is independent of the projection 
quantum numbers. This operator can be considered to act on a reduced basis, given by the 
complete set of basis vectors ||j(a)). In this form, the supcrblock wavefunction for target 
state j can be written 



II*) = fyUihaic*) llii(«i)> Ili2(a 2 )) , 

jlj20ClOC2 

over a product basis given by the Clebsch-Gordan series, 

||ji)®||i2> = || \ji-32\ )«•••© Nil +J2> ■ 
The reduced density matrix associated with this state is simply 



j 2 a 2 rj iai ;j 2 a 2 



(5) 



(6) 



(7) 



The matrix elements of the tensor product of operators T kl ® U k2 acting on different blocks 
are given by the Wigner 9j coefficients, 



3l 32 

k 1 fc 2 

fl 32 



3 
k 

f 



(8) 



(ilK)l|T fel ||jiK)}0' 2 K)ll^ fe2 lb2(« 2 )), 



where the [• • •] term is related to the 9j coefficient [12]. With this construction, all steps of 
the DMRG algorithm can be performed using only the reduced basis. The importance of this 
is that, unlike equation (1), there is no summation involved. The only essential difference 
from the standard DMRG formulation is the quantum number dependent 9j factor multi- 
plying each subspace. Thus, there is no significant computation penalty for using the SU (2) 
formulation, as long as the 9j coefficients can be calculated efficiently. In addition, for all two 
site interactions, the only two cases that appear are where one of the block operators in (8) 
is the identity operator, or when block operators are combined to form a rotational invariant. 
In both these cases, the 9j coefficient reduces to a single 6j coefficient. 

It is worth noting that in the SU (2) formulation, the basis vectors are exact eigenstates of 
total spin even after the truncation. This is not true, for example, if one attempts to force the 
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ground state to be in a particular total spin state by adding some suitably chosen multiple of 
S 2 to the Hamiltonian. Mixing of total spin states due to numerically near-degenerate states 
will still occur. Calculations involving long range interactions are also affected by the lack of 
explicit symmetries. Using a U(l) symmetric basis labelled by the z-component of spin only, 
interaction terms no longer transform as exact representations of SU(2) after a truncation. 
This can lead to situations where, even for a large number of kept states, the ground state is 
a broken symmetry Neel type state [13] and only converges slowly to an eigenstate of S 2 . It 
must be emphasized that this is purely an artifact of the DMRG algorithm when appropriate 
symmetries are not explicitly preserved. 

We now have a formulation of DMRG in which the states transform as 2j + 1 dimensional 
irreducible representations of SU(2). However, it is clear that the general formulation is es- 
sentially independent of the details of the SU (2) algebra - given an arbitrary compact global 
symmetry group the only modifications to the formulation is a different series expansion cor- 
responding to Eq. (6) and coupling coefficients from Eq. (8). For example, SO{4) is easy to 
utilize because one can make use of the isomorphism 50(4) ~ SU{2) x SU(2)/Z 2 , so that the 
6j and 9j coefficients are simply the product of two ST/ (2) coefficients. The component of 
the algorithm that is model dependent is rather small, consisting only of the reduced matrix 
elements of the single site operators, typically obtained via the Wigner-Eckart theorem. The 
Qj and 9j coupling coefficients of the algebra only appear in the construction of the blocks, in 
an identical way all models that admit the symmetry group. This separation of the physical 
aspects (the reduced matrix elements of the single site operators) and the geometric aspects 
(the coupling coefficients of the symmetry group) makes the method comparatively easy to 
apply to a range of models. This is the main advantage of the non-Abelian formulation over 
the IRF-DMRG. In the latter case, without a special effort to factorize the coupling coeffi- 
cients, the Boltzmann weights are rather complex quantities, especially for large symmetry 
groups. We have applied the non-Abelian DMRG successfully to various models with global 
symmetries SU{2), U{1) x SU(2) [16] and SO (A) [14]. The SU(3) case is in progress. 

The computational advantage of the non-Abelian construction is two fold: (1) each reduced 
basis element corresponds to 2j + 1 basis states of the old representation, thus the storage 
requirement for the block operators is reduced for an equivalent number of block states. (2) 
the superblock basis can be projected onto an exact subspace of arbitrary total spin. As 
well as reducing the size of the target Hilbert space, this greatly simplifies the calculation 
of excited states that have total spin less than the total spin of the ground state. This is 
very useful for investigating magnetic phase transitions [14]. For ferromagnetic target states 
(or more generally, target representations with a dimension greater than one), it is possible 
to calculate to first order the splitting of the degenerate states due to a symmetry breaking 
field, trivially in the case of a uniform magnetic field h (where the splitting is just hm, for 
m = —j, — j + 1, . . . , j), or in other cases by calculating the projection of the wavefunction 
and the symmetry breaking operator onto each z-component of spin using the Wigner-Eckart 
theorem (4). 

A model for which the non-Abelian formulation is eminently suited is the Hubbard model 
[17], 

H = 1 E ( c U c j> + H - c ") + U H (»*.T - I) ("U - I) ■ (9) 

<i,j>,a i ^ ' ^ ' 

The main feature of interest in the Hubbard model is the additional charge SU (2) pseudospin 

symmetry [18], generated by 1+ = E 1 (- 1 ) 1 4t c Ii' 7 ~ = £*( _1 ) ic »,i c i,T and F ~ = ^ 5( n *,T + 
rijj — 1). In the resulting reduced SO (4) basis, the Hubbard model contains only two basis 
states per site, a spinon of spin 1/2 and pseudospin zero, and a holon of spin zero and 
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pseudospin 1/2. The single site operators are 2x2 matrices over this basis. Table I shows a 
comparison of the ground state energy for the half-filled Hubbard model for a 60 site lattice 
with t = 1, U = 1, for the usual U(l) x U(l) basis of number of particles and ^-component 
of spin, the U(l) x SU{2) basis of number of particles and total spin S 2 — s(s + 1), and the 
5*0(4) basis of total pseudospin I 2 = + and total spin. For the case of half- filling, where 
the ground state is a spin and pseudospin singlet, the dimension of the representation D is 
equal to the number of basis states that would need to be kept to achieve the same accuracy 
using only U(l). Table I shows that the use of 5*0(4) symmetry gives an improvement of four 
orders of magnitude in the cumulative truncation error and the fractional error in the ground 
state energy, for virtually no increase in CPU time. The main contribution to the variance 
in the CPU times shown arises from differences in the number of matrix-vector multiplies 
being performed by the eigensolver, rather than any significant difference in the CPU time 
per matrix-vector multiply. 

Table II shows a calculation for a higher spin state, at half-filling with spin s = 5. In this 
case, the relative improvement from using SU{2) symmeteries is not as good. This has two 
causes. Firstly, the reduction in the dimension of the Hilbert space for total spin s — 5 versus 
z-component of spin s z = 5 is not as big as for the spin zero case. Secondly, s = is a special 
case in which the number of terms in the Clebsch-Gordan expansion (6) for the supcrblock 
is exactly one per block quantum number. For higher spin states, this is no longer true and 
the number of states in the superblock progressively increases as the target spin is increased, 
as each symmetry sector of the block basis appears multiple times in the superblock basis. 
Thus for a fixed number of states kept, the dimension of the superblock is much larger, with 
a corresponding increase in the CPU time. 

However we note that targetting a state of spin j is equivalent to inserting a non-interacting 
spin of magnitude j into the system and targetting a state of spin zero of the combined system 
+ non-interacting spin. Inserting this spin into the centre of the lattice, in between the system 
and environment blocks, is equivalent to targetting the spin j state directly. If however the 
non-interacting spin is placed at one end of the chain and integrated into one of the blocks (it 
doesn't matter which), the superblock dimensionality problem is avoided with only very minor 
loss of accuracy. This technique has been used before to target higher spin states in the IRF- 
DMRG algorithm [10]. The reflection symmetry is explicitly broken by the non-interacting 
spin, however the increased efficiency is well worth the loss of this symmetry. During the 
initial build sweep of the finite size algorithm, it is still necessary to target higher spin states 
because for a large enough non-interacting spin, there are not enough spins in the initial 
four-site block for the spin zero sector to contain any basis states [20]. For a target state of 
L sites and spin j, the actual target state during the build sweep when the lattice size is I 
sites (not counting the non-interacting spin), is j(l — l/L), rounded to the nearest permissible 
half-integer. Table III shows the improvement when this form of targetting is used. As this 
table shows, the additional superblock states that are included when the direct targetting 
method is used have negligible effect on the variational energy. 

We have extended the DMRG algorithm so that the block and superblock basis states 
transform as representations of an arbitrary compact non-Abelian global symmetry group, 
and demonstrated the improvement in accuracy for the Hubbard model utilizing spin and 
pseudospin symmetry. This is a true generalization of the conventional DMRG algorithm 
in that, if we instead use the coupling coefficients of U(l) instead of 51/(2) in Eq. (8), the 
original DMRG algorithm is recovered exactly. Thus optimizations such as efficiently storing 
the block operators [21], and transforming the obtained wavefunction to be the initial vector 
for the next DMRG iteration [22] apply to the non-Abelian case in a straight forward manner 
and were used in the current calculation. We have shown that, for the ground state of the 
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Table I - Comparison of (7(1) x 17(1), (7(1) x SU(2) and 50(4) basis for the groundstate of the 
half-filled Hubbard model for a 60 site lattice, at t = U = 1. Number of states kept m, dimension 
of the group representation D, energy E, fractional error in the energy, cumulative truncation error 
over the sweep 1 — a, CPU time in seconds per sweep. 



basis 


m 


D 


E 


(E-E g )/\E g \ 


1 -a 




CPU 


(7(1) x C7(l) 


100 


200 


-61.7484986435 


5.2xl0" 5 


5.3x10" 


■4 


10 


U(l)xU(l) 


200 


200 


-61.7514641444 


4.5xl0" 6 


4.8x10" 


■5 


41 


U(l) x 17(1) 


300 


300 


-61.7516910404 


7.9xl0" 7 


8.8x10" 


6 


110 


17(1) x 5(7(2) 


100 


226 


-61.7515581914 


2.9xl0" 6 


3.1x10" 


■5 


15 


(7(1) x 5(7(2) 


200 


468 


-61.7517319907 


1.3xl0" 7 


1.4x10" 


6 


64 


17(1) x 5(7(2) 


300 


716 


-61.7517389831 


1.4xl0" 8 


1.5x10" 


■7 


158 


50(4) 


100 


526 


-61.7517351742 


7.6xl0" 8 


8.4x10" 


■7 


18 


50(4) 


200 


1136 


-61.7517397636 


1.4xl0" 9 


1.5x10" 


■8 


71 


50(4) 


300 


1766 


-61.7517398448 


9.9xl0" n 


1.0x10" 


■9 


133 



half-filled Hubbard model, keeping only 300 SO{4) states is equivalent to keeping over 1700 
states of the U(l) x (7(1) basis of the original DMRG formulation. As the spin of the target 
state is increased, the accuracy improvement diminishes because the difference between the 
dimension of the Hilbert space of the total spin symmetry sector and the highest weight 
z— component of spin sector is reduced. Directly targetting a higher spin state is inefficient 
because the Clebsch-Gordan expansion implies that the dimension of the superblock will be 
much larger for the same number of block states. This inefficiency can be easily overcome by 
using a non-interacting spin to force the target state into the singlet symmetry sector. 

* * * 
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discussions on the IRF-DMRG algorithm. All calculations were performed on a desktop 
computer with a 500 MHz Athlon processor. 



Table II - Comparison of the (7(1) x (7(1), (7(1) x 5(7(2) and 50(4) basis for the lowest spin 5 
excited state for the half-filled Hubbard model on a 60 site lattice with t = U = 1. 



basis 




m 


E 


(E-E g )/\E g \ 


1 - a 




CPU 


t/(l) x 


(7(1) 


100 


-59.5701792131 


4.0xl0" b 


3.9x10" 


-4 


11 


U(l) x 


(7(1) 


200 


-59.5723270633 


3.6xl0" 6 


3.8x10" 


-6 


41 


(7(1) x 


(7(1) 


300 


-59.5725015232 


6.3xl0" 7 


6.8x10" 


-5 


102 


U(l) x 


5(7(2) 


100 


-59.5702795890 


3.8xl0" 5 


3.9x10" 


-4 


26 


C7(l) x 


517(2) 


200 


-59.5723402180 


3.3xl0" 6 


3.7x10" 


-5 


90 


U(l) x 


5(7(2) 


300 


-59.5725035338 


5.9xl0" 7 


6.8x10" 


-6 


207 


50(4) 




100 


-59.5723565660 


3.1xl0" 5 


3.5x10" 


-5 


32 


50(4) 




200 


-59.5725315037 


1.3xl0" 7 


1.4x10" 


-6 


127 


50(4) 




300 


-59.5725381508 


1.4xl0" 8 


1.5x10" 


-7 


297 
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Table III - Comparison of the U(l)xU(l), (7(1) x SU(2) and 5*0(4) basis for the lowest spin 5 excited 
state for the half-filled Hubbard model on a 60 site lattice with t — U =1, using a non-interacting 
spin to target the appropriate symmetry sector. 



basis 




m 


E 


(E-E g )/\E g \ 


1-a 




CPU 


(7(1) x 


SU{2) 


100 


-59.5702385716 


3.9xl0~ b 


3.8x10" 


-4 


14 


(7(1) x 


5(7(2) 


200 


-59.5723344203 


3.4xl0~ 6 


3.6x10" 


-6 


47 


U(l) x 


5(7(2) 


300 


-59.5725027479 


6.1xl0~ 7 


6.7x10" 


■6 


113 


S0(4) 




100 


-59.5723497975 


3.2xl0~ 6 


3.4x10" 


■5 


16 


SO(A) 




200 


-59.5725312667 


1.3xl0~ 7 


1.4x10" 


■6 


64 


SO(4) 




300 


-59.5725381253 


1.4xl0" 8 


2.2x10" 


-7 


148 
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